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1 Introduction 



The behavior of a quantum mechanical system in the semiclassical hmit strongly depends on the 
ergodic properties of the corresponding classical system. In particular, the eigenfunctions semi- 
classically reflect the phase space structure of the classical system and therefore they depend 
strongly on whether the classical system is chaotic or regular. In this work we are interested 
in the fluctuations of the wavef unctions, and in the correlations between the fluctuations in 
different regions which are induced by the classical phase space structures. In particular, we 
will consider the case of quantum billiards in a domain i7 C M^, which are described by the 
time independent Schrodinger equation 

(A + E)ilj{q) =0 for q G n\dn , (1) 

with Dirichlet boundary conditions, V'(q) = for q G dQ. For compact Q one obtains a discrete 
spectrum {En} of eigenvalues, < Ei < E2 < . . . , with associated eigenfunctions ipn £ -^^(^), 
which we assume to be normalized, i.e. \\ipn\\ '■= Jq |^n(q)P dg = 1. The corresponding classical 
billiard is given by the free motion of a point particle inside Q with elastic reflections at the 
boundary dfl. 

The amplitude distribution of an eigenfunction of a quantum mechanical system whose classical 
limit is chaotic is conjectured to become Gaussian in the semiclassical limit and numerical 
studies support this conjecture, see e.g. pj-Q. A more sensitive quantity is the local autocor- 
relation function [|l| which measures correlations between different points of an eigenfunction 

C^°"(x, 5x) := ^*(x - (5x/2)^/'(x + 5x/2) . (2) 

The crucial fact for the theoretical analysis of C^°'^(x, 5x), observed by Berry is that the 
autocorrelation function can be expressed as the Fourier transformation of the Wigner function 
(see eq. (J^) below) of ip, 

C'°=(x,5x)= / H^(p,x)e-'P^'' dp . (3) 



Hence information on the behavior of the Wigner function can be used to predict the behavior 
of the autocorrelation function, and since semiclassical limits of Wigner functions are concen- 
trated on invariant sets in phase space, see e.g. 0, it follows that in the semiclassical limit 
autocorrelation functions are determined by the classical phase space structure. For example, 
if the classical system is ergodic, the quantum ergodicity theorem |]6 HTl|] (roughly speaking) 
states that almost all quantum expectation values tend to the corresponding classical limit. 



One can show |T^ that for ergodic systems this is equivalent to the semiclassical eigenfunction 
hypothesis [|l3],0,[l^,|l5l, when restricted to a subsequence of density one. Using this result in 
(0) one gets Berry's result M that for chaotic billiards in two dimensions 



C'°^(x,5x)~--^Jo(v^|5x|) , (4) 

weakly as a function of x (for fixed 5x) as ^ 00, where E denotes the energy of the eigenstate 
•0 in (^. Equivalently we have 

hm Ci-(x, 6^/Ve) = -J— MM) . (5) 
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Numerical and experimental tests of this relation have been performed for several chaotic sys- 
tems 0-^,|T^ and show at finite energies noticeable fluctuations of the autocorrelation function 
around the high energy limit (H), especially for correlation distances larger than a few de Broglie 
wavelengths. These fluctuations have been studied further in [P^7|-pD|, where for small correla- 
tion length |(5x| a random model for the eigenfunctions of a chaotic system was used to predict 
the variance of these fluctuations, and for larger |5x| a formula involving closed orbits of the 



system has been derived. In |21|,|22| the path correlation function, which is an average of the 



local correlations along a given trajectory, has been introduced. The path correlation function 
is closely related to the autocorrelation function and for ergodic systems also tends asymptot- 
ically to a Bessel function (^. This path correlation function has been studied in for an 
hyperbolic octagon, and an expansion in terms of Legendre functions has been derived, which 
can be used to determine corrections to the leading Bessel part (H). 

The autocorrelation function in non-chaotic systems has attracted much less attention. The 
integrable case has already been discussed by Berry [Q, and the corresponding formula has been 
successfully tested for the circle billiard in [^. For a system with mixed classical phase space 



the autocorrelation function has been studied in [23|, in particular for irregular eigenfunctions 
an expansion of the Wigner function in polar coordinates has been used. 

In this work we are interested in the question how the universal limit (^) is reached, and how, in 
the case of mixed systems, further constraints on the classical motion are reflected in the auto- 
correlation function. For instance, if an eigenfunction is concentrated on an ergodic component, 



then by a generalization of the quantum ergodicity theorem [24], the Wigner function becomes 



equidistributed on that component, and this will determine the autocorrelation function. 
The paper is organized as follows. In section ^ we discuss some examples of the autocorrelation 
function for different eigenfunctions in systems with chaotic and mixed classical dynamics. In 
section |^ a general expansion of the autocorrelation function for eigenfunctions in billiards 
is derived, which allows a systematic study of their properties. It is an expansion in the 
correlation length \6x.\ which reflects the fact that the determination of correlations at larger 
distances needs classical information on finer length scales than for short range correlations. In 
section ^ it is shown that the correlation length expansion provides an efficient way to explain 
the fine structure of the autocorrelation functions of the systems studied in the first section. 
Of particular interest is that for chaotic systems deviations of the autocorrelation function 
from the quantum ergodic limit (|^) can be related to the rate of quantum ergodicity. In turn 
the autocorrelation function can be used to study the rate of quantum ergodicity on different 
classical length scales. 



2 Examples of autocorrelation functions 

For numerical computations as well as for theoretical considerations it is much more convenient 
to consider a smoothed version of the local autocorrelation function (H). Furthermore, as the 
eigenfunctions oscillate on a scale proportional to we rescale the autocorrelation function 

by this factor. Hence we will study the autocorrelation function in the form 

C,(x, 5x) := I p(x - q)^* (^q - ^) ^ (^q + 0Jj dq , (6) 
n 

where p is a positive function which determines the smoothing of the local autocorrelation 
function. In the literature (see the papers mentioned in the introduction) the mean is usually 
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taken over a small disk, which corresponds to taking the characteristic function of a disk for p 
in (^). However, nothing prevents one to consider the case p = 1, i.e. to take the mean value 
of the local autocorrelation function (H) over the whole position space. In terms of the Wigner 
function 



(27r)^ 



eW^*(q - q72)V^(q+ q72) dg' 



one obtains in this case 
C((5x) : = 



(5x \ / 5x 
l^(p, q)e-'P^"/^ dgdp = 



(p^|2g-ip5x/v^ 



dp . 



(7) 

(8) 
(9) 



This is a particularly good choice for the numerical computation of the autocorrelation function 
in billiards because it can be reduced to boundary integrals, see Appendix ^ The resulting 
formula reads 



C((5x) 



1 



8/E 



|q(sj -qi 



s') + ^x|yi(/E|q(s) -q(s') + 5x|)m*(s)m(s') dsds' , (10) 



where u{s) is the normal derivative of the normalized eigenfunction ip on the billiard boundary. 
This relation provides a very efficient method for the numerical computation of the autocorre- 
lation function. 

The systems for which we study the autocorrelation functions are the stadium billiard and two 
members of the family of limagon billiards, namely the cardioid billiard, and a billiard with 
mixed classical phase space. The stadium bilhard is proven to be strongly chaotic, i.e. it is 
ergodic, mixing and a i^-system [25, 26|. The height of the desymmetrized billiard is chosen to 
be 1, and a denotes the length of the upper horizontal line, for which we have a = 1.8 in the 
following. The family of limagon billiards is given by the simplest nontrivial conformal mapping 
of the unit circle ||2^,^ and can be parametrized in polar coordinates by p{ip) = 1 + e cos{(p) 
with ip G [— vr, vr] and e G [0, 1] denotes the family parameter. We consider the case e = 0.3 
which leads to a mixed dynamics in phase space. For e = 1 one obtains the cardioid billiard. 



which is also proven to be strongly chaotic pD|-pT|. The eigenvalues of the cardioid billiard 



have been provided by Prosen and Robnik and were calculated by means of the conformal 
mapping technique, see e.g. [^^. For the stadium billiard the eigenvalues and eigenfunctions 
have been computed using the boundary element method, see e.g. |^,^], and for the limagon 
billiard the eigenvalues have been computed using the conformal mapping technique and then 
the boundary element method has been used to compute the eigenfunctions (see |^ for details). 
For the high lying states in the limagon billiard the scaling method has been used |3^ . 
First we consider a "typical" eigenfunction in the cardioid billiard, see fig. In the plots we 
show 



Cir,e) = Cirei9)) , (11) 

where e(^) = (cos 6', sin 6'), as a function of r for three different values of 6. The quantum 
ergodicity theorem implies that there is a subsequence {rij} C N of density one such that 
Cnj{r,6) Joi'f') as rij — > oo with r fixed. This convergence is, however, not uniform in r. 



4 




Figure 1: Grey scale plot of |'?/'i277(q)P in the cardioid billiard with odd symmetry, 
where black corresponds to high intensity. To the right the autocorrelation function 
C{r,6), computed using ([T0|), is shown for three different directions 6 = 0,71/4 and 
tt/2. For comparison the asymptotic result C{r,6) = Jo{r) is shown as grey line. 




Figure 2: Autocorrelation function for the same state as fig. |l], but for a larger 
r- interval showing the non-universal behavior at larger r. The inset shows a mag- 
nification and the vertical bars indicate the places r = v^£^diam(f2, ^) from where 
on C(r, 9) = 0, due to the compactness of the billiard. 
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For the example shown in fig. ^ C (r, 6) fiuctuates, as expected for a "quantum ergodic" state, 
around the asymptotic result 

C(r,^)~ Jo(r) . (12) 

Actually, for an eigenstate with energy En we have C(r, 0) = for r > v'^diam(f2, 9), where 
diam(r2, 6') is the diameter of VL in the direction 6', as follows directly from the definition (^. 
This is illustrated in fig. |^ which clearly shows the non- universal behavior for larger r. 
In contrast to the case of quite uniformly distributed eigenf unctions one expects a stronger 
directional dependence of the autocorrelation function for localized eigenf unctions, such as 
scars [Q. One example is shown in fig. ^ where the eigenfunction shows localization along the 



shortest unstable periodic orbit in the cardioid. The corresponding autocorrelation function 
shows clear deviations from ([T2|). 

A class of eigenf unctions which show even stronger localization are the bouncing ball modes in 
billiards with two parallel walls, see e.g. [p9|,|^,^0|-|^. Fig. ^ shows for the stadium billiard an 



example of a bouncing ball mode, which localize on the so-called bouncing ball orbits having 
perpendicular refiections at the parallel walls and thus forming a one-parameter family. The 
simplest approximation is to consider them as a product of two sines, one in the x direction 
and the other in the y direction. In this case the autocorrelation function can be computed 
explicitly. For the odd-odd eigenfunctions 

'0n,,n^(a;, y) = l— sm{7m^x/l^) si\\{-nnyy/ly) (13) 



in a box B := [—Ix, Ix] x [—^y, ^y] one gets 

Cl:::nSr,e) = Fircosie)/VE,nx,lx)F{rsm{e)/VE,ny,ly) , (14) 




Figure 3: Grey scale plot of |'?/'„(q)p with n = 1817 in the cardioid billiard with odd 
symmetry. For the autocorrelation function C{r,6) one observes clear deviations 
from C{r,9) = Jo{r). 
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Figure 4: For the stadium billiard with odd-odd symmetry, a = 1.8, V'32o(ct) is a 
bouncing ball mode. The corresponding autocorrelation function is compared with 
the result C^^^{r,9), eq. (p!4D, obtained for a box, shown as dotted curves, which 
follow C{r, 6). Only for ^ = (full line) and 6 = n/A aX r ^ 17 are small deviations 
visible. 



where 

l-z/2 

I 



F{z,n,l) := xi~i,i]{z/2)j I sm{7Tn{x - z/2)/l)sm{nn{x + z/2)/l) dx (15) 



~l+z/2 

(l — — ) cos^nnz/l) + - — sm{Tmz/l) 
\ 21/ 2irn 



(16) 



and xi-i,i]{z) denotes the characteristic function of the interval [—/,/]. 

In fig. ^ we compare the autocorrelation function C{r, 9) for a bouncing ball mode in the stadium 
billiard with C\"^{r,9), eq. (H), and observe very good agreement. Mainly for ^ = some 



deviations are visible; these are understandable from the fact that in this case only correlations 
in the x-direction are measured, where the bouncing ball mode "leaks" outside the rectangular 
region. To take this into account one can determine an effective > 2a, by fitting sin^(7rx//^^) 
to 



i,rKx):= / \i^{x,y)\'dy . (17) 



For the case shown in fig. ^ this procedure leads to ~ 4 (whereas 2a = 3.6) and the 
corresponding autocorrelation function gives excellent agreement with the one for ip32o- 
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3 Expansion of the autocorrelation function 

In this section we derive an expansion of the autocorrelation function which will lead to an 
understanding of the directional dependence of the autocorrelation function observed in the 
last section. We start from the representation of the local autocorrelation function in terms of 
the Wigner function 

Cp(x,(5x)= /"/"p(x-q)Vr(p,q)e"'P^''/^dpdg . (18) 



Since the Wigner function is concentrated around the energy shell |p| = \/ and is furthermore 
even in p by time reversal symmetry, we get 

oo 27r 

Cp(x, 5x)= [ [ [ p(x - q)IV(p, q) dq e"'!^^! rdipdr + 0{E-^l^) 



p(x - q)IV(p, q) dg cos(|5x| cos((y9 - 9)) r dipdr + 0{E~^/'^) , 




n 



where we have used polar coordinates p = (|p| cosy^, |p| sin y^), (5x = (|(5x| cos6', |5x| sin6'). 
Because of the rescaling by y/E the factor q-^'p^^/^ is only slowly oscillating for p close to the 
energy shell, on which the Wigner function is concentrated. Therefore we get that the error is 
of order l/\fE, see Appendix ^ for a sketch of the derivation of this remainder estimate. If we 
now use that cos(rcosv9) is a generating function for Bessel functions 

oo 

cos(|(5x|cos(/?) = Jo(|5x|) + 2^(-l)'cos(2/(^)J2K|5x|) , (20) 

1=1 

we obtain 

oo 

Cp(x, 5x) = eo(x) Jo(|5x|) + 2 ^(-l)'6Kx, 9) MM) + 0{E-^'^) , (21) 



1=1 



with (setting r = |p|) 



oo 2iv 



6Kx,^):= j j j pi^-q)W{p,q) dqcosi2l{ip-9))rdifidr . (22) 



n 

The coefficients ^2i{'^,(^) can be further decomposed 

oo 2lT 



6/(x,0) = cos(2/^) j j y"p(x -q)Vr(p,q) dq cos{2lip) r dip dr 



n 

oo 2lT 



(23) 



+ sin(2/^) J J J P(x ~ q)W^(P; q) dg sin(2/v9) r dipdr 

Q 
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Recall that for an operator A with Weyl symbol A{p, q) the expectation value {ip, Aip) can be 
written as an integral over phase space of the symbol multiplied by the Wigner function of ip, 
see e.g. 



44 



(V^, AV^) = / / Wip, q) A(p, q) dp dq . (24) 



Therefore the coefficients in (^) can be interpreted as expectation values of certain operators 
A2/(x), -02/ (x) given as the Weyl quantizations of the functions 



^Ap, q) := p(x - q) cos(2/v?) 52z(p, q) = p(x - q) sin(2/9?) , (25) 



respectively, 



2n oo 



1 1 /^(P'^)'^(^"^)^°^(2M'^^'^^'^^^= ^^'^2«(x)V^) (26) 



n 

2tt oo 




W{p, q) p(x - q) sin(2/v5) r dr d(p dq = B2i{:s.)ip) . (27) 



n 



Note that the operators v42/(x) and -82; (x) depend on the parameter x. Since their symbols are 
smooth and homogeneous of degree zero in p they are classical pseudodifferential operators of 



order zero, see e.g. for the definition of pseudodifferential operators. So we finally obtain 



the following general expansion of the autocorrelation function 

C,(x,5x) = (^,lo(x)^)Jo(|5x|) 

00 

+ 2 l2/(x)^) cos(2/0) + B2ii^)ij) sm{2W)] J2ii\S^\) (28) 

1=1 
+ 0{E- 



1=1 

-l/2^ 



in terms of the expectation values of a sequence of bounded operators given as Weyl quantiza- 
tions of the symbols (^). Recall that the only approximation we have made was to insert for 
|p| in the exponent in eq. (|r^) the value at the energy shell y/E. 

Since the Bessel functions have the property that J2i(|(5x|) ^ for \6x.\ ^ 21, this representation 
is an efficient expansion for small |^x|; the larger |5x| becomes, the more terms of the sum have 
to be taken into account. Therefore it is desirable to have an estimate on the number of terms 
which have to be taken into account for large |5x|. The first, and largest, maximum of J2i{f) 
lies around r ~ 2Z, and close to it one has the expansion 

J2K2/ - zl^'^) = Ai(^) + 0(1/1) . (29) 

So the first peak becomes broader with a rate ~ l^^^ and therefore we have to take for large r 
approximately 
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terms in the sum over / into account; here z determines the error term. We refer to appendix 
y for a more detailed discussion. 



We would like to mention two works in which related results have been obtained. In p3[ a 
similar expansion was derived for the case that the eigenfunction is concentrated on an ergodic 
component of the phase space of a classically mixed system, however without extracting the 
Bessel function from the expectation values. For the case of a free particle on a surface of 
constant negative curvature an expansion of the path correlation function in terms of Legendre 
function was derived in P]. In the special case of averaging over the whole billiard (i.e. p = 1) 
the path correlation function should be for ergodic systems the same as the autocorrelation 
function. 

The correlation length expansion has various possible applications; some of them will be 
discussed and illustrated in the next section. In particular, the expansion leads to a prediction 
for the asymptotic limit of the autocorrelation function in different situations. More precisely, 
consider a subsequence of eigenfunctions {"^Anj jjeN for which the corresponding sequence of 
Wigner functions converges weakly to a measure v on phase space. Such a measure v is called 
a quantum limit, and it is an invariant measure of the classical flow. 

If a sequence of eigenfunctions {V^nj jjeN converges to a quantum limit, the correlation length 
expansion for the autocorrelation function (p8|) shows that the corresponding sequence of au- 
tocorrelation functions converges as well and their limit is obtained by substituting in (|28|) the 
expectation values of A2/(x) and B^ii^ by their corresponding classical limit. Explicitly, this 
gives 



Cj°^^*(x,(5x) = AoJod^xl) + 25^(-l)'[A2Kx)cos(2/^) + S2Kx)sin(2/^)] J2/(|5x|) , (31) 
where 



A:= j Adu . (32) 

As we will discuss in section for ergodic systems the terms A21 and B21 vanish for i^^ — >• 00, 
and with = 1 we recover (p!2D. 



4 Applications of the correlation length expansion 
4.1 Direct comparison 

In the numerical examples we have studied the autocorrelation function in the case p = 1, 
which allows for an exact computation of the autocorrelation function using the representation 
(PUP, which is much more efficient than a direct computation of the autocorrelation function 
by its definition, eq. (|]). In this case the general expansion (|28|) gives the representation 

00 

C(r, 9) = Mr) + 271 ^^(-1)' [a2i cos(2/^) + b2i sm{2W)] J2i{r) + 0{E-'/^) , (33) 
1=1 

where the coefficients 02/ and b2i are the Fourier coefficients 

a2i = — / /(v?) cos{2l(p) dip h')i = — I I{(p) sm{2l(p) dip , (34) 

TT J n J 
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Figure 5: Comparison of the autocorrelation function C{r, 6) for tpigoj in the stadium 
bilhard (full curve) with the expansion (|33|). In particular for small r the agreement 
is excellent, whereas for larger r small differences become visible. 

of the radially integrated momentum density [^,^|, 

oo 

:= j |^(re(^))|Vdr , (35) 



where e(y9) = (cosy?, sin </)). Also for I{(p) a representation in terms of a double integral of the 
normal derivative function is available [^. Taking the symmetries into account, one can show 



that for the odd eigenfunctions in the limagon billiards and the odd-odd eigenfunctions in the 

stadium billiard all b2i vanish, so only the cosine terms remain in ( P5| ) and (PSI). 

First we will test the influence of the error term 0{E~^^'^) in eq. (^) for computations at finite 



energies. To that end we use the exact quantum I{(p) in eq. (|5^. In fig. ^ the autocorrelation 
function C{r,6) for four different angles 6 is compared to (0). In particular for r not too 
big the agreement is excellent. Only for larger r do small deviations become visible, which go 
to zero for higher energies and r fixed. One should remark that for any r > the effective 
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Figure 6: Angular dependence of the autocorrelation function C{r,6) for different 
r. Shown are the results for ^/'eooo in the cardioid billiard with odd symmetry. The 
full line is the result for C{r,6) using (0), the dashed line shows the result of the 
expansion (|33|) , the full grey line is the value of Jo{t) and the dotted horizontal lines 
show the variance Jo{r) ± A^/^, see eq. (p7|). 
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integration region in eq. (H) is reduced by the factor 



vol (n n l](r//E, 6) 



where VL{r/yE,6) is the set Vt shifted by the vector r/yE (cos ^, sin ^). Incorporating this 
factor leads to an improvement in the agreement of the expansion with the exact autocorrelation 
function at larger r. 

Instead of looking at the dependence of the autocorrelation function C{r,9) for fixed 6 and 
varying r, it is also interesting to keep r fixed and consider the angular dependence. For a 
"chaotic" eigenfunction in the cardioid billiard some examples are shown in fig. ^ The result of 
the expansion ( |5BD is in good agreement with the exact result. For larger r the autocorrelation 
function C{r,6) oscillates more strongly around Jo{r). For even larger r we observe clear 
deviations of the expansion from the exact result (not shown). For comparison the variance 
of the autocorrelation function around the prediction Joir) for a random wave model |T^ in 
leading order 

AV.= r^V''i (37) 



is shown and good agreement is found. 
4.2 Localized eigenfunctions 

For a state strongly localized on an periodic orbit of length we have (either in the semiclassical 
limit, or as a crude model at finite energies) 



/(<^)~^5Z^%-¥'*) , (38) 
where /^^ are the lengths of the segments of the orbit with direction <^j. Thus we get 

asz = ^ ^7, cos(2/(/)i) = ^ ^7, sin(2/v9i) , (39) 

which therefore using (^) gives a prediction for C(5x) for such states, namely 

C(5x)~l5^/,,cos(|5x| cos(6' - (pi)) . (40) 



Notice that in the presence of symmetries all symmetry related directions have to be taken into 
account in eq. (p8|). For this simple model one can determine the autocorrelation function more 
directly by using (H) 

C(5x) = jj W^(p,q)e^P^^dpdg = j |t^(p)|VP^'' dp 

(41) 



2n 



I{(p) cos(|5x| cos(6' - (f)) dip + 0{E 



-1/2^ 
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inserting (^) directly gives (^). 

In fig. 1^ we compare the limiting behavior ( |1D| ) with the autocorrelation function of a high 
lying eigenstate in the limagon billiard. The state localizes on the (stable) orbit of triangular 
shape. Up to r ^ 10 the agreement is very good; for larger r the autocorrelation function 
of the eigenstate shows deviations from the asymptotic behavior. Notice that the state has a 
much higher energy than the other examples. At lower energies the agreement is not as good, 
because the region in phase space on which the state localizes is broader. This in turn implies 
that its corresponding radially integrated momentum distribution I{(p) also has broad peaks, 
which is not accounted for properly by the ansatz (|38D. However, when considering states of 
this type with increasing energies, a clear trend to the asymptotic result (|40|) is observed. 
This simple model has also been tested for a scarred state in the cardioid. However, the 
agreement is limited to a qualitative description for up to r 2. This is understandable 
in view of the observation (see [^, fig. 8a)]) that for a scarred state the radially integrated 
momentum distribution I{(p) shows quite large fluctuations, and also in the considered case the 
direction ip = tt /2 is not clearly pronounced. As these fluctuations essentially correspond to 





Figure 7: High lying eigenfunction {E = 367984.82 . . . , approx. 47788**^ eigenfunc- 
tion of odd symmetry) in the limagon billiard {e = 0.3), which localizes on the stable 
orbit of triangular shape. The autocorrelation function for three different directions 
is compared with the 5-model, eq. (|4D|), shown as dashed line using the directions of 
the stable orbit. 
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the random "background" fluctuations of tlie state, a simple ansatz to model this behavior is 
C{r,9) = {l-a)Jo{r) + aY^l^,cos{\6^\cos{9-^i)) . (42) 

It turns out that one can vary a such that quite good agreement of this model with the exact 
autocorrelation function is obtained, see fig. |] where a = 0.22 (for all directions). Depending 
on the direction 6 the "optimal" value for a does vary, which already indicates the limitations 
of this simple model. To get a better agreement a more precise description of I{(p) for scarred 
states is necessary. In particular, this should also lead to an understanding of the energy 
dependence of a which is expected to go to zero in the semiclassical limit. Notice, that the 
structure of the autocorrelation function is quite similar to the one for ipisn shown in fig. ||. 
Another case, for which we obtain much better agreement, is for an eigenfunction localized on 
an invariant torus. In such a case the expectation values, eqs. (pGf), (|27|) , tend to the mean of 



the classical observable over the torus, see eqs. (pTD, ([321). Fig. |^a) shows for the limagon billiard 
the eigenfunction and the corresponding Husimi Poincare representation ^,|4^; see for 



a more detailed discussion and the formula which has beens used. Also shown in the Husimi 
plot are the points of some orbits. Using an initial condition on the torus we can determine 
the classical angular distribution /classical ^^^^ j^^ ^^lis has a singularity due to the caustic of the 
torus we show in fig. ^) a binned distribution together with the corresponding quantum radially 




Figure 8: For a scarred state ("^7147 of odd symmetry) in the cardioid billiard the 
autocorrelation function is compared with the simple model (|42|) for a = 0.22. 



15 






Figure 9: Grey scale plot of V'sose for the limagon billiard with e = 0.3 together with 
the corresponding Husimi plot, for which in addition some orbits are shown. In c) the 
radially integrated momentum distribution /sosel^') and the corresponding classical 
distribution /classical ^^-j ^^^q torus are shown. In d)-f) the exact autocorrelation 
function is compared with the expansion of the autocorrelation function, eq. (^), 
using /classical j^j^ -J different angles 6. 
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integrated momentum distribution /soselv')- There is qualitative agreement between these two 
curves in the sense that smoothing /classical ^^-j (describes the mean behavior of the quantum 
-^3056 (v^)- Of course, the classical distribution cannot describe the (quantum) oscillations visible 
for /3056(v^)- It turns out, see fig. ^)-f), that already this simple model leads to surprisingly 
good agreement between the exact autocorrelation function and the expansion (|33|) computed 

using /classical -J _ 



4.3 Autocorrelation function of irregular states in mixed systems 

In classical systems with mixed phase space regions with regular and regions with stochastic 
behavior coexist. It is conjectured ||50| that correspondingly the quantum mechanical eigenfunc- 
tions split into regular and irregular ones, respectively, living semiclassically on the correponding 
parts of phase space. This has been confirmed numerically for several systems, see e.g. [pTH55[|. 
Consider now a sequence of eigenfunction ip^^ which localize on some open ergodic domain D 
in a system with mixed phase space, then almost all the expectation values (ipn^Aipn^) tend to 

the mean of the corresponding classical observable A over this domain D [^]. Therefore 
using (|33[) we get in the limit E ^ oo for the autocorrelation function of such a sequence 

oo 

C^^'\^, 5x) = A^"" Jo(|5x|) + 2 ^(-1)' [ A^''(x) cos(2/e) + 5^''(x) sin(2/^^)] M\6^\) . (43) 

1=1 

Instead of computing A2i^ and B2i^ directly, we can also use a typical trajectory of the ergodic 
component to determine the corresponding classical /classical j^^-j 

/'="(^)= limy5^/.5(^-^.) , (44) 

where / is the total length of the trajectory and (fi is the direction of the i-th segment having 
length /j. Then we use (|33| ) to get a prediction for the autocorrelation function. 
However, we observe that even quite high lying states do not yet localize on the whole chaotic 
component. Instead they are confined to smaller subregions due to partial barriers in phase 
space. Fig. [TO|a) shows an example of a high lying state in the limagon billiard {e = 0.3) In 
fig. p!0[b) the corresponding Husimi function is plotted, which clearly shows the localization on 
a chaotic subdomain (the whole irregular region is much larger). If D is a open region in phase 
space, then the corresponding classical distribution of the momentum directions is given by 

/^"(V') = I Xd{p{^), q) dg , (45) 

where p{^p) = (cosy?, sin yj). One can show that in terms of the projection V of D on the 
Poincare section this equation can be reduced to 

Jl{s,p)5{(p - 0(s,p)) dsdj9 
'""""'^^^ = }l(„),s,p <^«) 

V 

= , (47) 

J l{s,p) ds dp ' 

V 
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Figure 10: Autocorrelation function for a high lying irregular state {E = 
1002754.70 . . . , approx. 130516**^ eigenfunction of odd symmetry) in the limagon 
billiard with e = 0.3. In b) the Husimi representation on the boundary is shown 
together with an approximate boundary (full line) of the region V on which the 
state localizes. The resulting classical momentum distribution /classical j^^^-j -g gj^own 
in c) as full line and compared with the radially integrated momentum distribution 
J^™(y9) of the state in a) and a smoothing of this, /qni'Smoothedj^^^^^ shown as dashed 
line. In d)-f) the autocorrelation function C{r,6) of the eigenfunction is compared 
for three different directions with result of the expansion (^3p using /classical 
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where l{s,p) is the length of the orbit segment starting in the point {s,p) G V with direction 
4>{s,p) and in the second equation p{s,ip) = p{(p)t{s), with t(s) denoting the unit tangent 
vector to dD in the point s. Furthermore dQ'{(f) := {s G dQ \ p{ip)n{s) < 0}, where n{s) 
denotes the outer normal vector to dD in the point s, is the subset of dfl where the vector 
p{if) points inwards. For the numerical computation we have used p7|) because we just have 
to deal with a one- dimensional integral to compute the (p dependence, and also compared to 
( |46| ) no binning of /classical j^^^-j -g ngcggga^-y. 

After these general remarks on the computation of /classical (jggcj-jbg j^qw we compute 
the relevant quantities to determine the autocorrelation function for the state shown in fig. p!0|a). 
To describe the projection T> of the domain D in phase space, we use an approximation of the 
boundary of P by a splines, which are shown in the fig. ^b) as full curves. Then we use 
eq. (0) to determine the corresponding /classical ^^^^ shown in c) as full curve. Of course the 
radially integrated momentum distribution I^"^{ip) of the eigenstate shows strong fluctuations, 
but the smoothing /qm.smoothedj^j^-j -g ^^jj described by /classical j^^^-j^ although the agreement 
is not perfect. Using /classical j^^^-j employ the expansion (^31) to get a prediction for the 
autocorrelation function for states localizing on V, which is compared in figs. p!0|d-f) with the 
exact autocorrelation function. Up to r 10 we get quite good agreement, whereas for larger 
r deviations become more visible. This shows, that the eigenfunction has more structure than 
accounted for by /classical ^^^^ jg enough in the semiclassical limit. 

For higher energies the states tend to localize on the full ergodic region, and then /classical 
can simply be computed using (|^) by averaging of a typical trajectory in D. One should 
emphasize that the agreement has to be compared with the agreement of the autocorrelation 
function for ergodic systems with ([12|) as the prediction eq. ( ^31) only takes into account the 
classical limit. This has been studied in (in the case of averaging the local autocorrelation 
function over a small disk), where in particular for [^, Fig. 13b] very good agreement has been 
found. 



4.4 Ergodic systems and the rate of quantum ergodicity 

If the classical system is ergodic and ipn- is a quantum ergodic sequence of eigenfunctions, then 
for j oo 

{^jjnj , A2I (X) ^jjnj ) ~ = (5;o (48) 

(7/'..,S2/(x)Vn,) ~ fi^ = . (49) 



Thus using the expansion (pSD we again get (|l^) for — 00. Deviations from this universal 
behavior are then determined by the rate of quantum ergodicity. In order see this it will be 
convenient to remove the angular dependence by taking the mean over all angular-directions in 
C(5x). Since by eq. 



2n 



2n 

we consider the second moment. 



^ / C„(r, d) dd = Mr) + 0{E-'/') , (50) 
ivr J 



2n 



i- j[C^{r,e)-Ur)fde , (51) 



19 



where Cn{r,6) denotes the autocorrelation function of ipn- Inserting the expansion (|33|) of 
C„(5x) leads to 

oo 

^l{r) = 27r'Y,{al^ + blJ[Mr)]^l + 0{E-'/')) . (52) 
1=1 

In fig. |TT] we compare a^{r) for an eigenfunction in the stadium billiard with the expansion 
(|52|) . For small r we get excellent agreement and some deviations become visible in the plot for 
r > 20. The inset shows a plot of the difference up to r = 20. It is surprising that even though 
for large r the amplitudes do not match anymore, still the expansion gives the right oscillatory 
structure. 

If we take the mean of (|52D over all eigenfunctions up to energy E, we get 

^(E,r) := J2 -nir) (53) 

OO 

= 2^' E ME) E (4,n + 4,n) [Mr)r (1 + 0{E-'/')) (54) 

1=1 ^ ' E„<E 

Remarkably, together with eqs. (^8]) and (|52|) this shows that the rate of quantum ergodicity 
can be studied in terms of the autocorrelation function. Particularly interesting is that the 




Figure 11: Comparison of the second moment crlgQ^i^r) of the autocorrelation func- 
tion, eq. (^Tj), with the expansion (^) for the stadium billiard. The inset shows the 
difference for r G [0, 20]. 
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observables in the expansion (|2^) become more and more oscillatory with increasing Z, so by 
varying |5x| one can determine the rate of quantum ergodicity on different length scales. 
A prediction for the behavior of ct^(£', r) follows from |56|, where it is shown that (under suitable 
conditions on the system) in the mean 

,55) 

for any pseudodifferential operator A of order zero with symbol A. Here A denotes the mean 
value of A, and ali{A)/VT is the variance of the fluctuations of 

T 

i J A(p(t),q(t)) dt (56) 



around A. So if we insert (ESl) into (|5l) we obtain 



9 CO 



^'iE,r) ~ +^i(i?20] IMr)?^ . (57) 

A detailed study of the rate of quantum ergodicity in terms of the autocorrelation function, 
i.e. via eq. (|53D, and a comparison with the semiclassical expectation ( [57| ) will be given in a 
separate paper. 



5 Summary 

We have discussed the autocorrelation function for eigenstates of quantum mechanical systems, 
and its relation to the behavior of the classical system. For billiards we have derived a formula 
for the autocorrelation function of an eigenfunction in terms of the normal derivative on the 
boundary (p!OD, which enables an efficient numerical computation. 

Our main result is the correlation length expansion of the autocorrelation function (^8]), which 
provides an efficient expansion for small correlation length, where only a small number of 
terms enter. Moreover, it provides a tool for understanding the behavior of the autocorrelation 
function for different types of eigenfunctions in terms of their semiclassical limit. 



The coefficients in the correlation length expansion (28) can be computed in terms of the 
radially integrated momentum density. Even though it is based on an approximation, our 
numerical study shows very good agreement with the corresponding exact results; only for 
large correlation lengths deviations become visible. As the expansion coefficients have to be 
determined just once for a given eigenfunction, this is also a numerically efficient method to 
compute the autocorrelation function. Similar expansions can be derived in higher dimension 
and for more general systems (e.g. systems with potential and magnetic field), but then the 
Bessel functions have to be modified in order to reflect the structure of the energy shell of the 
classical system. 

We applied the expansion of the autocorrelation function to different types of eigenfunctions, 
and showed that it provides a good tool for the understanding of their autocorrelation func- 
tions. In systems with mixed phase space regular states concentrated on tori and irregular 
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states have been successfully treated. For chaotic system the fluctuations of the autocorrela- 
tion functions around the leading term are shown to be connected with the rate of quantum 
ergodicity. Moreover by varying the correlation length the autocorrelation function is shown to 
be an interesting new tool to measure the rate of quantum ergodicity on different length scales. 
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Appendix 

A Autocorrelation function in terms of normal deriva- 
tives on the boundary 

We will give a derivation of the formula ( [T0| ) which provides an expression of the autocorrelation 
function C{Sx.) in terms of the normal derivative. Let ip{ci) be a solution of the Helmholtz 
equation with Dirichlet boundary condition on dQ, 

(A + k^)ilj{q) = , ^(q) = for q G , (58) 

where we have defined k = \/E, and let 

u{s) := n(s)V^/'(q(s)) (59) 

be the outer normal derivative of ip on dQ, where s parameterizes dfl in arclength. It is well 
known that 



j/>'oWq-qWI)..Wd.= (*<l) 

4 J for q e i7 



(60) 



and furthermore 



j Mk\q-q{s)\)u{s)ds = . (61) 



an 



Let p{t) be a smooth cut-off function with 

fl for t< 2diam(fi) 
p(kt) = { - ^ V (62) 

^ ^ [0 for t>3diam(n) ^ ^ 

where diam(f2) denotes the diameter of Q. Then we have for q in some neighborhood of Q 

^icL) = -\ j p{k\ci-q{s)\)Yo{k\ci-q{s)\)u{s)ds (63) 

an 
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and obtain 



C(5x) = y"^*(q)V'(q + 5x) dg = jj Kp{Sx, s, s')u* {s)u{s') ds ds' , (64) 



anxdn 

with 



Kp{5^,s,s') = ^ J p(fc|q-q(s)|)yo(fc|q-q(s)|)yo(fc|q-q(s') + (5x|)dg 
= ^ y"p(^|q|)>o(A:|q|)ro(fc|q + q(s)-q(s') + 5x|)dg . 



(65) 



Due to the factor p(A;|q — q(s)|) this integral is absolutely convergent. We now use Grafs 
addition theorem [^] 

y„(^|q + Aq|) = |S«^'«^<ll)^'(*=l''l'™^('^' for|q|<|Aq| 

\E,aliWq|)^i(*|Aq|)cos(/4,) for|q|>|Aq|, 

where Aq = q(s) — q(s') + (5x and ip is the angle between Aq and q. Introducing polar 
coordinates in the integral in (|65|) and using (^) gives 

|Aq| 

K,{6^,s,s')=^ J Yo{kr)Mkr)r drFo(A;|Aq|) 

' oo (67) 

+ p(A;r)ro(A;r)Fo(/i;r)r dr Jo(A;|Aq|) , 

|Aq| 

where we have furthermore used that p{kr) = 1 for r < |Aq| by (|62D. The first integral is 

|Aq| ^ 

Yo{kr)Mkr)r dr = [yo(A;| Aq|) Jo(A;| Aq|) + ri(A;| Aq|) Ji(fc| Aq|)] , (68) 







see e.g. |H3|, and for the second one partial integration gives 



p{kr)Y„{kr)Y„{kr)r dr = -^!L [y„(fc| Aq|)y„(*-| Aq|) + y,(fc|Aq|)n(*:|Aq|)] 

(69) 

k f 

-- p{kr)[Yo{kr)Yo{kr) + Yi{kr)Yi{kr)]r^ dr . 

|Aq| 

Note that since p' has compact support the second integral is over a finite interval, and for 
s,s' e dQ, 5x G the lower limit of the integral, |Aq|, is outside the support of p', hence the 
second term on the right hand side of eq. (pO) is constant. So we get 



Kp((5x, s, s') = K{5x, s, s') + Rp{5x, s, s') (70) 
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with 

K{6^, s, s') = [Fi(A;| Aq|) Ji(fc| Aq|)Fo(fc| Aq|) - Y,{k\Aq\)Y,{k\Aq\) Mk\ Aq\)] (71) 

and 

Rp{6^,s,s')=CMk\Aq\) (72) 
with C constant and by (|6l|) this term gives no contribution to C(5x). Using a Wronsky 



determinant of Bessel functions we can simphfy K{S:s., s, s') further 



ir(5x, s, s') = !^^ri(A;| Aq|) [Ji(fc| Aq|)ro(A:|Aq|) - ri(A:| Aq|) Jo(fc| Aq|)] 
which gives the final result. 

B Remainder estimate 

In this appendix we sketch the derivation of the remainder estimate in equation ([l9|) . We start 
by representing the integral as an expectation value, see (p^), 

oo 2tt 

j j y"p(x-q)W^(p,q) dg'e^"l''"l™^('^-^)/^rd^dr = (7/',v47/') (74) 
n 

where A is the Weyl quantization of the symbol 

a(p, q) := p(x - q)e'^l'"l '"'^''-'^ . (75) 
The basic idea is to find a decomposition of the operator A, 

A = Aq + (v^ - /E)Ai + R (76) 

where Aq has Weyl symbol 

ao(p,q)=p(x-q)e^l^''l=°^(^-^) (77) 

and the remainder R satisfies 

\\R\\< CE-^'^ . (78) 
If we assume the decomposition (|76|) and take the expectation value of both sides, one gets 

(7A,A7/>) = (7/>,AoV') + (^,/?V') , (79) 
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where {y/—A — VE)^p = has been used. In terms of the symbols eq. (|79D is the desired result, 
see m). 



oo 2n 



Q 



j j y Vr(p, q)p(x - q) dg'e'''l'^"l^°^(^-^)/^r dv^dr 

oo 2lT 



(80) 



n 



Let us now show, that the decomposition ([76|) is basically a quantization of the Taylor expansion 
of the symbol a(p, q) around |p| = V^, 



a(p,q) = ao(p, q) + (IpI - v^)ai(p,q) . (81) 

Quantizing this classical decomposition yields ([75| ) with R given as Weyl quantization of 

^(P,q) = (IpI - v^)ai(p,q) - (IpI - v^)#ai(p,q) (82) 

since the Weyl symbol of {^/—/S. — \/E)Ai is (|p| — V^)#a,i(p, q) with # denoting the symbol 
product, see e.g. Since is a constant we have 



^(p, q) = |p|ai(p,q) - |p|#ai(p, q) 



(83) 



and this is a function which is bounded and of order 0{E ^/^), and all its derivatives are 
bounded and of order 0{E^^/'^), too. So by the Calderon Vallaincourt theorem the estimate 
(|7S|) follows. 

C Estimating the Bessel sum 

In this appendix we determine how many terms in the sum ( PB] ) have to be taken into account 



such that the remainder is smaller than some given error 6. From (^6]) and (pT]) it follows that 
for fixed if) 



\{ij,A2i{^)i,)cos{2ie) + (V^,52Kx)^)sin(2/^^)| < C . 
Thus if we split the sum 

oo 

5^(-l)'[(^,l2/(x)^) cos(2/^) + (^,52Kx)V') sin(2/^)] J2i{\5^\) 



(84) 



1=1 



m—l 



(85) 



^(-1)'[(V^, A2i{^)ij) cos(2/e) + B2i{^)i:) sm{2W)] J2i{\S^\) + i?m(|^x|) 



1=1 



we get for the remainder 



\RUr)\<Cj2\Mr)\ . 



(86) 



l=m 
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Therefore we have to estimate the sum over Bessel functions 

oo 

T.\Mr)\, (87) 

l=m 

and determine its dependence on m and r. The asymptotics in the transition region 

M2l-z{2iy/')^j^M{2'/'z) , (88) 
gives that J2i{r) is monotonically increasing for r < 21, such that for r < 2m 

1 — m 7 — m ^ ' 



(89) 



Defining z by 



r = 2m - zm^/^ , (90) 



we obtain 



^^±j^. ( 2{l-m) /m^ 

oo , / ^, / \ 1/3 



~ ^ (Z + m)V3 + m)V3 + ^ [j^ ) ) ^^^> 



E^Ai(J, + .)+0(™-V=) 

1=0 ^ ^ 



where we have furthermore used that for large m only the terms with / -C m contribute, because 
the Airy function is exponentially decreasing for positive arguments. The Euler McLaurin 
formula then gives 



oo °° 

E ( + =) = / ^ Ai ( + .) d/ + 0(m-'/») 

/— n 





oo 



(92) 



Ai(a;) dx + 0{m-^'^) , 



And so finally we arrive at 

„ oo 
oo „ 

E \Mr)\ = 2 / ^^(^) + 0(m-i/3^ . (93) 

The function Ai(a;) dx is monotonically decreasing, so for a given (5 > we can define a 

z{5) by 

oo 

^ j Ai{x) dx = 5 , (94) 
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Figure 12: For the bounds 6 = 10~^ and 6 = 10^^ of the sum over Bessel functions 
(P^) the result of the exact computation of m(r, 6) (full curves) and the asymptotic 
result (^) are compared. The asymptotic results approaches the exact one slowly 
from below with a rate 0{r~^^^). 



and then (^) defines together with (0) a function m(r, 6) such that 



|J2/(r)| = (5 + 0(r-i/3) . (95) 

l=[m(r,5)+l] 



By solving (|90D for large r, we see that we have to take approximately 

m(r,5) -- + -(-) (96) 

terms in the sum (|28| ) over Z into account such that the error is 5 + 0{r~^l^\ 

For instance, if we require 5 = lO""^, then (5^) gives ^(5) = 4.359 . . . ; for 5 = 10~^ one gets 



z{S) = 7.925 .... In fig. |T2| we show for these choices of z the asymptotic result (|96D compared 
to the exact computation, corresponding to (P7|). The asymptotic results approaches the exact 
one slowly from below; in the plotted region a constant offset by two compared to (|96D gives a 
good bound for m(r, 6). 
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